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ABSTRACT 

We show that accretion disks, both in the subcritical and supercritical accretion rate 
regime, may exhibit significant ampHtude luminosity oscillations. The luminosity time 
behavior has been obtained by performing a set of time-dependent 2D SPH simulations 
of accretion disks with different values of a and accretion rate. In this study, to avoid 
any influence of the initial disk configuration, we produced the disks injecting matter 
from an outer edge far from the central object. The period of oscillations is 2 — 50 
s respectively for the two cases, and the variation amplitude of the disc luminosity 
is 10"^* - 10'^^ erg/s. An explanation of this luminosity behavior is proposed in terms 
of limit cycle instability: the disk oscillates between a radiation pressure dominated 
configuration (with a high luminosity value) and a gas pressure dominated one (with a 
low luminosity value) . The origin of this instability is the difference between the heat 
produced by viscosity and the energy emitted as radiation from the disk surface (the 
well-known thermal instability mechanism) . We support this hypothesis showing that 
the limit cycle behavior produces a sequence of collapsing and refilling states of the 
innermost disk region. 
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1 INTRODUCTION 

This work continues our studies on the occurrence of the 
Shakura and Sunyaev instability dShakura fc SunvaevllQT^ 
in the a-disks when the radiation pressure dominates, i.e. in 
the so-called A zone. 

The problem of the existence and outcome of the Shakura- 
Sunyaev instability is important in accretion disc physics 
because it affects the models and their time behavior. 
In general, the outcome of the Shakura-Sunyaev instability 
is guessed to be the formation of a hot cloud around the 
internal disc region, in which comptonization could happen 
(Shapiro ct al. 197^. 

A critical point of this scenario is the typical timescale re- 
quired by the disk to leave the collapsed 'dead' state. Time 
dependent analytical models of the disk evolution in this 
post collapsed phase are very difficult. Numerical simula- 
tions are therefore important and essential tools to obtain 
some indication of the outcome of this evolution. 
Recently some authors have investigated this problem 
through the numerical appr oach. Szuszkiewicz and Miller 
JSzuszkiewicz fc MilleJ Il99'<t) found that a slim accretion 
disc model with low viscosity {a = 0.001) and a luminosity 
higher than O.OSLb shows a thermal instability which gives 
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rise to a shock- like structure near to the sonic point, leading 
to the disc disruption. They found no limit-cycle behavior, 
probably, according to their own conclusions, because of the 
not str ong enough advection. The same Szuszkiewicz and 
Miller ISzusz kiewicz fc Millei|[l998ft also performed numeri- 
cal simulations of accretion disc models with high viscosity 
(q = 0.1) and obtained a limit-cycle behavior. They sim- 
ulated the disc evolution for several cycles and, for M = 
0.06Af_B and a central object of 10 solar masses, found a pe- 
riod of the cycle of about 780 s. In both papers they reported 
the results concerning a vertically integrated disc model, 
with no consideration of acceleration in the vertical direc- 
tion. The same authors dSzuszkiewicz fc Miller! l200jl per- 
formed finally a numerical study of an accretion disc model 
at high viscosity (q = 0.1) with a vertically integrated treat- 
ment of acceleration in the vertical direction and a diffusive 
form for the viscosity instead of the aP prescription used 
in their previous works. Also with this more refined model 
they found a limi t-cycle behavior. 

Nayakshin et al. jNavakshin et al^ll985^ used a limit-cycle 
model to explain the luminosity variability of the micro- 
quasar GRS 1915-1-105. Their model is different from the 
one used by Szuszkiewicz and Miller. The essential differ- 
ence regards the viscosity prescription. Szuszkiewicz and 
Miller used the standard Shakura-Sunyaev one or the more 
refined (but fundamentally equivalent) diffusive formulation. 
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In their ID simulations the discs osciUate between two sta- 
ble states, one at high luminosity and the other with a 
much lower emission. These two states are the standard 
gas-pressure dominated one and the radiation pressure dom- 
inated one (note that this last state is stable in the slim ac- 
cretion disc model). If the accretion rate difference between 
the high and low states is very large, as in GRS 1915-1-105, 
the high state should last a very short time. But in reality the 
GRS 1915-1-105 has a high state lasting for a long time, even 
more than the low state. Therefore the limit-cycle model in 
slim accretion discs cannot explain the time behavior of this 
source. So Nayakshin et al. used a particular viscosity law 
which produces a high stable state of larger duration. With 
this model, they explained the gross observational features 
of GRS 19154-105. 

Janiuk et al. jjaniuk et aljl20o3l also tried to describe the 
GRS 1915-1-105 behavior in terms of a limit-cycle model. 
They adopted the standard a-viscosity prescription, but in- 
cluded in the model the effect of a corona surrounding the 
disc and a vertical outflow. With half of the energy dissi- 
pated in the disc, they obtained outbursts whose amplitude 
and duration are consistent with the GRS 1915H-105 data. 
Teresi et al. jTeresi et al.1 1200^^ have clearly shown that, 
at intermediate accretion rates, accretion discs with A zone 
suffer a collapse but after a rather long time they show a 
flaring activity with an intervening refllling phase of the A 
zone. 

We point out that our simulations differ from the 
Szuszkiewicz and Miller ones since we produce real 2D disks 
with true vertical motion. No 'ad hoc' prescription is re- 
quired to include physical vertical effects. The new degree 
of freedom given by the Z motion has many consequences, 
which we discuss in section 4. 

Furthermore we point out that all previous simula- 
tions start from a full disk existing at time t — 0. A 
typical drawback of simulations involving large disk sec- 
tors is the uncertainty in the initial model structure. In- 
deed the analytical models lack a reli a ble vertical str ucture 
jBisnovatvi-Koean fc Blinnikovl IT977I: iHubenvl ll990D . The 
disk luminosity and other disk parameters can oscillate for 
long time before reaching a steady configuration and it can 
be difficult to discern between real instability oscillations 
and transient spurious oscillations, that could last a long 
time. 

The simulations we show here differ from our previous 
ones and from many other authors since the disc is gener- 
ated "ab initio". We inject matter at a large distance from 
the central object. The injected gas has low temperature and 
keplerian angular momentum. Its evolution is due to the ac- 
tion of the viscous stress, whose a value is given. In this way 
the disk evolves smoothly through a series of equilibrium 
states, avoiding the problem of the transient spurious oscil- 
lations and the influence of the initial configuration on the 
final result. It is clear that also for these simulations a long 
integration time (of the order of the viscous drift time) is 
required. The numerical Smoothed Particles Hydrodynam- 
ics technique allows to integrate up to so long times. Let 
us note that -in general- a lagrangean code, as SPH is, is 
better suited to capture convective motions than eulerian 
codes. With the same spatial accuracy (cell size equal to 
particle size) the SPH particle motion is tracked with great 



accuracy, i.e. the particle size may be large but its trajectory 
can be still determined 'exactly'. 

Our results suggest that the Shakura-Sunyaev model 
can be used to explain the luminosity variability shown by 
many sources. The aim of this work, however, is not to find 
an explanation of some sources behavior, but simply to see 
what happens to the standard disc structure during its time 
evolution. 

The paper is structured as follows: in section 2 we 
remind the physical model, in section 3 we describe the 
adopted numerical method, in section 4 we report on the 
simulated cases and the obtained results, presenting and 
commenting some figures, in section 5 the physical aspects of 
the simulations results are discussed and in section 6 we ex- 
pose conclusions and astrophysical implications of our work. 



2 THE PHYSICAL MODEL 

The time dependent equations describing the physics of ac- 
cretion disks are well known. We used the lagrangean form 
of them in a cylindrical reference system and in the approxi- 
mation of local thermal equilibrium (LTE) between gas and 
radiation ( .Mihalas fc Kleiniil98^'l . They include: 
mass conservation 



Dp ,. ^ 

— — — ~p div V 
Dt ^ 

radial momentum conservation 

P-p^ = + pgr + (dlV a)r + fr 



vertical momentum equation 



Dt 
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energy equation 
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where aij , a , is the viscosity stress tensor and 9 is the 
body force per unit mass (acceleration), 
angular momentum equation 

= -2n-^ + -—,,p—) + - — — {r\p—) 5 
Dt r p oz az pr-^ or or 

Here Q, is the local angular velocity, is the comoving 

derivative, Erad is the radiation energy per unit volume, / 
is the radiation force per unit volume, given by: 



/= P F 



F= 



E is the radiation flux, given by: 

V Erad 



3p{k + gt] 



(6) 



(7) 



k and cjt are the free-free absorption and Thomson scat- 
tering coefficients, given by: 
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k = CkpT 2 [cm g 



(8) 



(Kramers' formula, valid since our temperatures are 
well beyond lO" K) 

with Cfc = 2.26 • 10^*, and 

ar = OA[cm^g'^] (9) 

A is the angular momentum per unit mass, E = e+ -^^^ 
is the total internal energy per unit mass, including gas and 
radiation terms. 



The components of a that we have considered in our 
calculations (since they are the ones that play an important 
role in accretion discs) are the r-(j) one, given by: 



(Jr0 = vpr 



dr 



and the 0-2 one, given by: 
d{Q.r) 



dz 



(10) 



(11) 



v = aVaH is the kinematic viscosity, a is the viscos- 
ity parameter of the Shakura-Sunyaev model, Va is the local 
sound speed, H is the disc vertical thickness and the other 
terms have the usual gas dynamic meaning. 
The gravitational field produced by the black hole is 
given by the well - know n pseudo-newtonian formula by 
iPaczvnski fc Wiital Jl980l) : 



Eqrav — 



{R 



GM R 
R9Y R 



(12) 



where R is the position vector of the point in which the 
field is evaluated, with a modulus given hy R = Vr^ + z^, 
Rg is the Schwartzschild gravitational radius of the black 
hole, given by: 

2GM 



Rg — 



(13) 



and M is the black hole mass. 



We adopt the local thermal equilibrium approxima- 
tion for the radiation transfer treatment. However this 
assumption does not affect our conclusions. 



domain and the same procedure to search the near neighbors 
of each particle. Therefore applying the usual procedure for 
the evaluation of any smooth function in the point {ri,Zi) 
we have: 

/(r0^//(O^.(.-r')|^.r'. 



E/0 



(14) 



where = {rk,Zk)- So for the density we have the simple 
expression that identically satisfies the continuity equation 
in the cylindrical form: 



N 



(15) 



Rewriting the fundamental equatio ns in the formul ation 
more suitable for the SPH evaluation iMonaiha^liii), and 
applying the previous criteria we have the following expres- 
sions; for the radial (r) momentum we obtain: 



\ Dt J, 



r ^ r, \p} p2 



dm, 

dri 



(16) 



where 11 is the artificial viscosity pressure. 



The vertical (z) momentum satisfies: 



( — ) 

\ Dt J, 



dzi 



(17) 



For the energy equation we based our implementation 
on the following procedure. Let us call U = Ejrad j^^j^ 1^2 ^ 

-P*°' = Prad + Pgaa ,and remind that F= - ^ V {Prad + Pgaa) 

is the total force per unit mass due to gas and radiation, then 
the first three terms of the energy formula can be trivially 
put in to the SPH form alism according to standard prescrip- 
tions iMonaghanlll992fl . So we obtain: 



dt ^ \ 



V ■ g + V ■ F 



" /ptot ptot\ ^ 



3 THE NUMERICAL METHOD 

We set up a new version of the Smoothed Particles Hy- 
drodynamics (SPH) code in cylindrical coordinates, for axis 
symmetric problems. We remind that SPH is a lagrangean 
interpolating method. Recently it has been shown it is equiv- 
alent to finite eleme nts with spa rse grid nodes moving along 
the fluid flow lines (iDilts 111996'). For a detailed account of 
the SPH algorithm see iMon aghan (1985). For cylindrical 
coordinates implementation s ee iMolteni et all J|l998)) and 
IChakrabarti fc Moltenil (|l99rf ). Our code includes viscosity 
and radiation treatment. 

The basic point for our cylindrical geometry approach is 
simply to assume a usual kernel function but depending di- 
rectly on the radial (r) and vertical (z) variables, and there- 
fore retaining the usual normalization factor and width. Now 
pseudoparticles are small tori of mass druk = 2iv g^rkdrkdzk 
. In this way we may use the same Cartesian grid in the (r, z) 



+ 



here Vik=Vi — Vk- The fourth term 



(18) 



can be 



written, using the same method of SPH evaluation for the 
V V term of the continuity equation, as: 



N 

E 

fc=i 



'Pj±Pk.' 

K 2 , 



Sik -ViWik 



(19) 



where we have symmetrized the density term, and where 

S=v ■ CT, Sik=Si — Sk', with this procedure the SPH en- 
ergy equation conserves exactly the total energy. 

The total thermal internal energy ^^^"^ + e is recovered 
by subtraction of the kinetic energy and then the ratio be- 
tween .5^^ and e is given requiring the LTE condition. 

In cylindrical coordinates the particles masses nii must 
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be replaced by ^ and a further Vk term appears in the 
terms coming out of the divergence expressions, so we have 
for example: 



E 



rftot 

p 



— \ — — 

' Pi 



\ ^ mfe 



— tcyi 



Pt 



Pi 



— >cy! 



where Vik = i^iVvi — r^Vr^) r + (u^. — t)^^) z , and r and z 
are the radial and Z versors. 

To derive all previous expressions we neglect the contri- 
butions to the integrals from the boundary of the integration 
domain. The artificial viscosity pressure flij is formulated as 



pi] 



(20) 



with the averaged quantities 

_ -|- Cj 



Pij = 



pi] — 



+ 



pi + P] 



1% = (n -Tjf + {Z^ - Zjf 



T) = O.l/l 



a and fi are the artificial viscosity coefficients used to damp 
out oscillations in shock transitions, c here denotes the sound 
speed. 

Since our aim is to simulate accretion discs it was es- 
sential a correct treatment of the tangential velocity and 
its diffusion due to the viscosity. We integrate explicitly the 
viscosity diffusion term; the cylindrical SPH-version of the 
dif fusion term of equ ation @ is given following the criteria 



bv lBrooksha^ lll994) : 



QiQ] 



D ^ 

^i] 



where 



(21) 



(22) 



The formulae for cylindrical geometry are similar to the 
cartesian ones and the most relevant changes are : (i) the 
mass of a particle appears divided by its distance from the 
z-axis (ii) mutual velocity difference Vj — Vi between two 
particles must be replaced by the more sophisticated term 
{riVi - rjVj)/ri. 

The force Fr differs from Fr while F^ - - = F^ .. for 
particles at the same radial coordinate. This difference in the 
force is due to the geometry. Angular momentum is exactly 
conserved. The statement druk = 2n QkVkdrkdzk is needed 
only for the derivation of the formulae and the particles in 
the simulations may have the same mass or not; obviously 
the density is no more directly proportional to the number 
of particles per unit area as for the case of particles having 
the same mass. 

To integrate the energy equation we adopted the split- 
ting procedure. In the LTE condition the radiation energy 



density changes according to the well known diffusion equa- 
tion given by: 



-div F = V • f V Eraa 



dt \3pntot 

where ktot = k + ar- 

In cylindrical coordinates r, z: 



(23) 



dErad C d 



dErad\ ^cd_ ^ r dErad I j.^^^ 



dt r dr \ Sputot dr J r dz \ Spntot dz 
where c is the light speed 



The SPIf- version of the radiation transfer term is given 
following the criteria by Brookshaw ( 1994) . The cylindrical 
coordinate version is given by: 



1 m, f Ei - Ej^ \ ^ Kii 



( — \ ^ ± \ " ^ 

\ dt ) i Ti ^ Tj 



Pi 



^]J^-^^W,, (25) 



where for clarity we did not put the subscript rad in 
Erad and where: 



+ 



■ipif^toti SpjKtot 



-z,){26) 



This formula can be obtained by the same procedure 
explained by Brookshaw, but taking into account that -in 
cylindrical coordinates- the particles masses are defined as 
rrik = 2tt pk Vk Arfc Azk, that explains the further division 
by rj in the term ^ . 

W e used a variable h procedure iNelson fc Papaloizoul 
Il994fl . In our procedure, in order to have a not too small 
particle size (and therefore not too great CPU integration 
times), we put a fioor for the h values: h is chosen as the 
maximum between the value given by the variable h proce- 
dure itself and 1/10 of the disc vertical half thickness. So we 
have nearly 10 particles along the disc half thickness even in 
the collapsed region. 

The boundary conditions of the simulations are not 
fixed, though we produce an infiow at a certain radius, gen- 
erating new particles with fixed density and temperature 
every time a circular zone around the injector position be- 
comes empty. As the SPH particles move around, the simu- 
lation region follows the form assumed by the disc and the 
values of the physical variables at the boundary of the disc 
are the values that characterize the boundary particles at a 
certain time. 

The spatial extension of the initial configuration is decided 
by establishing a radial range of physical interest and a ver- 
tical extension given by the disc thickness of the Shakura- 
Sunyaev model. The values are given in the next section. 
For radiation, the boundary conditions we used are based 
on the assumption of the black-body e mission and partic - 
ularly on the Brookshaw approximation iBrookshawlll99.tf) . 
At every time step boundary particles are identified by geo- 
metrical criteria (the particle having the maximum absolute 
value of z in a vertical strip of radial width given by /i is a 
boundary particle). The boundary particle loses its thermal 
energy according to the formula given by Brookshaw (that is 
an approximation of the diffusion equation at the single par- 
ticle level), that states the particle cooling rate proportional 
to iQT)/h^, where Q = [AacT^) / {ipntot)- 
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In all our simulations the boundary particles never 
reach an optical thickness lower than 10. 

4 THE SIMULATIONS PERFORMED 

We performed several simulations, the ones commented here 
had the following parameter values: 

a) a = 0.1, M = 0.15, domain Ri - R2 = 3 - 100, 
h = 0.25; 

b) a = 0.1, M = 2, domain iii - i?2 = 3 - 200, h = 0.5; 

where M is in units of Me and Me is the critical ac- 
cretion rate. For all cases the central black hole mass is 
M = 10 Mq. The initial spatial resolution we adopt is h. 
The reference units we use are Rg for length values, Rg/c 
for time values and Ltheor = 0.06Mc^, the theoretical lu- 
minosity for an accretion disc around a not rotating black 
hole, for the luminosity values. In case 'a' Ltheor = 1.54 • 
10^^ erg sec~^; in case 'b' Ltheor = 2.05 • 10^® erg sec~^. 
We have chosen these units because the simulations results 
obtained with certain values of the parameters M, M and 
a, if given in terms of adimensional units, can be easily gen- 
eralized to other systems. 

The radial domain of each simulation has been chosen with 
the aim of including in the simulation a sufficiently wide por- 
tion of the radiation dominated zone, the so-called A zone 
(in the case 'a' the entire A zone is included in the simula- 
tion). The h values above reported are the initial ones. They 
have been chosen in order to have a good spatial resolution 
at the injection radius. The variable h procedure guarantees 
then an equally good resolution in the inner disc regions. 

For case 'a' we stress that our results have been obtained 
simulating a full disc including A, radiation dominated, and 
B, gas pressure dominated, zones. The presence of the B 
zone, that is theoretically stable and that we see stable in our 
simulations, guarantees in general the numerical accuracy 
of our study and allows to clearly identify the A zone as 
responsible of the oscillations. 

In this section we want to show the changes that occur 
in the main properties and physical quantities of the disc due 
to the instability and the consequent limit-cycle behavior. 
When the instability arises the disc undergoes a collapse 
phase, with a strong lowering of its vertical thickness. Fig. 
1 makes evident the effect of this phenomenon, showing, for 
the case 'a', the disc configuration reached at the end of the 
collapse phase, characterized by a very small Z height in 
the innermost region (r < 13Rg). Note that the Z scale is 
graphically amplified in the figure. 

In this state the mass accretion rate is no more uni- 
form throughout the whole disc. In fact, in the innermost, 
collapsed disc zone M has a value lower than the one be- 
fore the collapse, whereas at the outer boundary of this zone 
it assumes the value of the outer, not collapsed region, i.e. 
the unperturbed value. So mass is forced to enter the col- 
lapsed zone at a rate larger than the one at which mass falls 
into the black hole. Due to this accretion rates difference, 
the collapsed disc zone is refilled and consequently reaches a 
configuration of much larger vertical thickness (comparable 
to the one of the unperturbed state). In fig. 2 we show the r- 
z profile of this new structure, together with the boundaries 
of the corresponding ShakurarSunyaev disc, determined by 
calculating the disc vertical thickness (at all the r values in- 
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Figure 1. The r-z profile of tfie disc of case 'a' in its low state is 
shown at the time t = 0.52155 10^ Rg/c ■ Every SPH particle is 
represented by a small cross. On the x axis the r values in units 
of Rg are represented. On the y axis the z values in units of Rg 
are represented. 
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Figure 2. The r-z profile of the same disc case 'a' in its high 

state at the time t = 0.52233 10^ Rg/c is shown. The solid line 
represents the equivalent Shakura-Sunyaev model. 

side the simulation radial range) from the Shakura-Sunyaev 
one-dimensional model with the same accretion rate and a. 
This comparison of the 2D simulated model and the ID the- 
oretical one is shown in order to make evident the agreement 
we obtained between the results of the simulations and the 
canonical disc model. We will discuss more deeply this point 
in the section 5. 

Moreover, the difference between the two states (col- 
lapsed and refilled) is not only in the value of the disc ver- 
tical thickness. In the unstable disc region the temperature 
of the collapsed state is lower than the one of the refilled 
configuration. As a consequence of that, the ratio between 
radiation and gas pressure is changed from a state to the 
other one: though in the unstable region the disc always re- 
mains radiation pressure dominated, during the collapsed 
state the ratio Prad/Pgas is much lower (close to 1) than in 
the refilled disc. In fig. 3 we show the comparison between 
the radial profiles of the ratio Prad/Pgas in the two states, 
which makes evident the great lowering of this ratio in the 
transition to the collapsed state within the unstable disc re- 
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Figure 3. The ratio Prad/Pgas at the times t = 0.52155 10^ and 
t = 0.52233 lO'^ is shown. On the X axis the r values in units of 
Rg arc represented. The configuration at the carher time exhibits, 
in the collapsed zone, a much smaller ratio Prad/Pgaa than the 
configuration at the later time. 



5.19 5.2 5.21 5.22 5.23 5.24 5.25 5.26 5.27 5.28 5.29 
time 

Figure 5. The time behavior of the disc luminosity is shown. On 
the X axis the time values in units of Rg/c are represented. On 
the y axis the luminosity values in units of L^hnor = 0-06Mc^ = 
1.54 • 10^* erg sec~^ are represented. 



Shakura-Sunyaev moael 
High state 
Low state 




Figure 4. The temperature radial profile at the times t = 
0.52155 lO'^ and t = 0.52233 lO'^ is shown. On the x axis the 
r values in units of Rg are represented. On the y axis the tem- 
perature values in Kelvin are represented. The configuration at 
the earlier time exhibits, in the collapsed zone, a much smaller 
temperature than the configuration at the later time. The solid 
line represents the temperature radial profile for the equivalent 
Shakura-Sunyaev model. 



gion, and in fig. 4 a similar comparison for the temperature 

profiles is shown, making clear that the unstable region is 
cooler in the collapsed state than in the refilled one. In fig. 4 
is also shown the temperature radial profile (in the disc mid- 
plane) of the corresponding Shakura-Sunyaev model. Here 
also we show the good agreement of our simulations with 
the canonical disc model. 

In these figures, and in all the figures of this paper in 
which physical quantities are plotted versus r, we show the 
values regarding all the SPH particles: to each particle (of 
radial coordinate, say, r) corresponds in the figure the point 
(r, Q), where Q is the value of the physical quantity that 
we are plotting calculated in the particle position. 
The temperature difference also produces different luminosi- 
ties associated to the two configurations. So the disc lumi- 



nosity oscillates between the two states and we can observe 
the limit-cycle behavior typical of the thermal instability. 
In fig. 5 we show the time variation of the disc luminosity, 
from which the oscillatory behavior is clear. In this figure 
only a time window of the luminosity variation regarding 
the whole history of the disc is represented. The time units 
are, as said, Rg/c. To obtain the time values in seconds it is 
sufficient to multiply the values in the figure by 10"'* (the 
value of Rg/c in seconds for a black hole of 10 solar masses). 

What can be noticed, in particular, from this figure is 
the shape of the time variation curve. A single oscillation 
starts with the disc luminosity L that increases very steeply; 
then, after having reached a maximum, L decreases more 
slowly (with an exponential-like behavior) until a value close 
to the initial one is reached. 

We also evaluated the gas velocity field, finding a significant 
difference between the radial speeds ( Vr) in the two states: 
in the unstable region the refilled disc has a higher radial 
speed (with a large spread) than the collapsed one. This is 
what we can expect considering that the refilled disc is more 
luminous (since it is hotter) and therefore the accretion rate 
of its inner region is larger with respect to the collapsed disc. 
A larger accretion rate can be the effect of a larger radial 
speed. In fig. 6 the radial profiles of Vr in the two states are 
shown. 

It is evident from this figure what we have said above 
and also that in the refilled state the radial speed is often 
positive, besides very high. If a large radial speed is present 
in both infiow and outflow directions, as clearly shown in 
the figure, we can argue that there is a gas circulation and 
not only a large net accretion rate. 

We highlight that this result excludes the possibility that, in 
2D, a Shakura-Sunyaev disc can show a regular radial flow. 
This simplifled one-dimensional picture is destroyed by the 
presence of convective and circulatory motions in the r — z 
plane. These motions can be considered as a new turbulence, 
different from the one that gives rise to the a- viscosity: the 
equations of the disc dynamics contain a viscosity term (the 
a-term) that is physically considered the result of a sup- 
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Figure 6. The gas radial speed is shown versus r in the two 
states (the high and the low one) of the disc. On the x axis the 
r values in units of Rg are represented. On the y axis the radial 
speed values in units of c are represented. 
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Figure 7. The r-z profile of the disc of case 'b' in its low state is 
shown at the time t = 0.1016946 lO'^ Rg/c . Every SPH particle 
is represented by a small cross. On the x axis the r values in units 
of Rg are represented. On the y axis the z values in units of Rg 
are represented. 



posed turbulent motion, but a disc simulated according to 
these equations develops another turbulent motion (that can 
be supposed to give rise to an additional viscosity). Our re- 
sults about circulation and convection and our hypothesis 
that these 2D motions can be at the origin of an additional 
viscosity agree with the analytical study by Kippenhahn and 
Thomas on con vective and circulatory flows in thin radiative 
accretion discs jKippenhahn fc ThomaJl982l) . The large ra- 
dial speed has also to be considered the reason for which 
the thermal instability causes the disc collapse and not its 
expansion. The local perturbative approach in itself allows 
to conclude that a small temperature deviation from the 
equilibrium state, an increase as well as a decrease, grows 
exponentially in time. Therefore the result should be, with 
the same probability, an expansion or a collapse. What we 
observe, instead, is that collapse is strongly preferred: in 
each cycle, initially the inner zone reduces largely its verti- 
cal thickness, then it swells reaching a thickness value not 
much larger than the equilibrium one. Our hypothesis is that 
what lacks in the local perturbative approach is the radial 
drift of matter, and therefore energy, due to the advective 
motion. This radial flow, carrying away thermal energy from 
a disc region at a certain r before the expansion instability 
has developed at that radius, inhibits the local thermal en- 
ergy growing and therefore the disc expansion. 

For the case 'b' we show the r-z proflles of the disc 
particles in the two states in flgs. 7 and 8. In particular, 
in fig. 8 we also show the boundaries of the corresponding 
Shakura-Sunyaev disc, from which the reader will be able to 
see the good agreement, in the vertical thickness, between 
our simulation and the canonical disc model. 

The reader will notice that the whole disc is geometri- 
cally thinner in a state than in the other one. The unstable 
region is no more a small zone close to the black hole, as 
in case 'a', but extends throughout the entire simulation 
radial range. This is due to the higher (supercritical) accre- 
tion rate, that makes the radiation pressure dominated zone 
much wider than in case 'a'. The extension of the unstable 
region is also apparent from the comparison between the ra- 
dial profiles of the ratio Prad/Pgas in the two states, shown 
in fig. 9, and from the similar comparison for the tempera- 
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Figure 8. The r-z profile of the same disc case 'b' in its high 
state at the time t = 0.1167006 10^ Rg/c is shown. The solid line 
represents the equivalent Shakura-Sunyaev model. 

ture profiles, shown in fig. 10, where the temperature radial 
profile (in the disc midplane) of the corresponding Shakura- 
Sunyaev model is also included (in order to show also here 
the good degree of agreement between simulations and ID 
models) . 

In fact, in both figures the two profiles associated to 
the two states are significantly different in the whole disc: 
approximately up to 200-Rg, where the disc mass is injected, 
the geometrically thinner conflguration is cooler and less 
radiation pressure dominated than the other one (the disc 
remains however radiation pressure dominated) . So what we 
said about the 'a' case is also valid in the 'b' one: we can 
conclude that there is a limit-cycle oscillation between two 
different disc states, one thinner, cooler and therefore less 
luminous and the other one thicker, hotter and therefore 
more luminous. The luminosity oscillation is shown in fig. 
11. In this figure the luminosity time behavior is represented 
during the entire formation and evolution of the disc from 
the time t — 0, when particles begin to be injected in a 
totally empty space. 
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Figure 9. The ratio Prad/Pgas at the times t = 0.1016946 lO'^ 
and t = 0.1167006 lO'^ is shown. On the X axis the r values in 
units of Rg are represented. The configuration at the earlier time 
exhibits a smaller ratio Pradl Pgas than the configuration at the 
later time. 



Figure 11. The time behavior of the disc luminosity is shown. 
On the X axis the time values in units of Rg/c are represented. On 
the y axis the luminosity values in units of Ltf^^^r = 0.06Mc^ = 
2.05 ■ 10^^ erg sec~^ are represented. 



Shakura-Sunyaev model - 
High state 
Low state 




Figure 10. The temperature radial profile at the times t = 
0.1016946 lO'^ and t = 0.1167006 lO'^ is shown. On the x axis the 
r values in units of Rg are represented. On the y axis the temper- 
ature values in kelvin are represented. The configuration at the 
earlier time exhibits a smaller temperature than the configuration 
at the later time. The solid line represents the temperature radial 
profile for the equivalent Shakura-Sunyaev model. 



in the equation above is the same for the whole disc region, 
i.e. the energy density E{r, i) is uniform in the considered 
region. 

It is easy to notice that the oscillation shape in the 
'b' case is more similar to th e light curves obtained in ID 
simulations (see, for example, 'Watarai & Mineshige' ('20031), 
[Navakshin ot al. (1985) and Szuszkiowicz & Miller (2001)) 
than the luminosity behavior of the 'a' case. Our hypothesis 
to explain this fact is that in the 'b' case the unstable zone is 
much wider. So the mass to be unloaded by the disc to pass 
to the collapsed state is larger. Because of this the refilled 
state duration is greater and the oscillation shape shows 
a luminosity maximum characterized by a width approxi- 
mately equal to the minimum duration. In fact, we have, 
in the 'b' case, two states with luminosities different by a 
factor 7 whose durations are about 25 a each one. 
In the 'a' case, instead, when the luminosity maximum is 
reached the fight curve starts immediately its descending 
(exponential) phase: the disc holds for a very short time 
the refilled state, probably because the disc gets soon rid of 
the small amount of matter contained in the small unstable 



It is evident from this figure the rather different shape 
of the time variation curve with respect to the analogous 
curve for the 'a' case. Here the luminosity L increases ap- 
proximately as steeply as it then decreases. So there is no 
exponential-like behavior as there is in the 'a' case. We guess 
that this difference in the luminosity behavior is due to the 
much larger extension of the unstable zone in the 'b' case. 
In fact, the cooling of a given disc portion is probably gov- 
erned by a law closer to the exponential one if the energy 
density of the considered disc region is fundamentally uni- 
form throughout the region itself. If we indicate with E(r, t) 
the local energy density at the radius r and the time t in the 
disc, the cooling law at the position r is = —kE{r, t), 

that has for solution the exponential form for E(r,t). 
It is obvious that this argument holds for the total energy of 
an entire disc region only if the variable E{r,t) to consider 



Finally, we studied the radial behavior of the Mach 
number M = — . 

We present the results of this analysis in the figs. 12 
and 13, regarding respectively the collapsed and the refilled 
state. For clarity the figures show only the innermost disc 
region. It is clear from these figures that the disc has a sonic 
point, positioned nearly at r = IO-R9 in the collapsed state 
and (very approximately) at r = 15i?g in the refilled one. 
From the external boundary to the sonic point the radial 
flow is subsonic, whereas from the sonic point to the inter- 
nal boundary we have a supersonic flow. Though here our 
data are strongly scattered, we can say that our 2D simu- 
lated models reveal a transonic region at radii larger than 
in ID calculations. For a comparison with ID models on 
this aspect of the disc dynam ics, one can see, for example, 
ISzuszkiewicz fc Milleil 1)1998}) , where the sonic point is given 
around r = SRn- 
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Figure 12. The radial profile of the Mach number M = for 
the disc of case 'b' in its collapsed state. Every SPH particle is 
represented by a small cross. On the x axis the r values in units 
of Rg are represented. On the y axis the Mach number values are 
represented. 
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Figure 13. The radial profile of the Mach number M = ^ for 
the disc of case 'b' in its refilled state. Every SPH particle is 
represented by a small cross. On the x axis the r values in units 
of Rg are represented. On the y axis the Mach number values are 
represented. 



5 DISCUSSION 

In this section we discuss three items: 

1) confirmation by a true 2D approach of the limit-cycle 
behavior produced by the thermal instability; 

2) differences between the results obtained by the two 
approaches ID and 2D; 

3) theoretical considerations about the main features of 
the luminosity time behavior; oscillations amplitude, typical 
times. 

1) As we said in the first section, the limit-cycle 
behavior due to the Shakura-Sunyaev instability has al- 
ready been shown as result of ID time-dependent simula- 
tions iSzus zkiewicz & Miller 199 7. 1998, 2001; Janiuk ct aL, 
|2002|) . In this work we confirm the existence of the limit- 
cycle behavior using a 2D approach, obviously closer to the 
physical reality of accretion discs. Moreover, the 2D ap- 
proach allows to reveal aspects of the accretion flow that 




Figure 14. The gas velocity field in the radial range between 
SSiJj and 53-Rg is shown for case 'b'. The disc state is the high one. 
On the X axis the r values in units of Rg are represented. On the 
y axis the z values in the same units are represented. The arrows 
represent the velocity vectors, with their lengths proportional to 
the speed values. 



Velocity field of the particles in a portion of the disc 




Figure 15. The gas velocity field in the radial range between 
33iJg and 53i?g is shown for case 'b'. The disc state is the low one. 
On the X axis the r values in units of Rg are represented. On the 
y axis the z values in the same units are represented. The arrows 
represent the velocity vectors, with their lengths proportional to 
the speed values. 



cannot be simulated by the ID methodology. Convection and 
circulation are the main ones. Also, convective motions are 
supp osed to reduce the thermal instability ( Shakura ct all 
Il97al . Just in presence of a signiflcant convective flow, that 
our 2D simulations reveal, we confirm the existence of the 
thermal instability. 

2) In figs. 14 and 15 we show the gas velocity field in 
a given disc region for the two disc states. The considered 
case is the 'b' one. 

These figures make clear the different feature of the 
two states (high and low) as regards the convective motions 
in the disc. From fig. 14 we can argue the relevant pres- 
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ence of convective and circulatory motions in the high state, 
whereas from fig. 15 it is evident the poor convection in the 
low state. 

These phenomena are important because they affect the 
time values that characterize the luminosity behavior. We 
note that the light curve we obtain is different from the 
curves usually obtained by ID simulations. The luminosity 
behavior we find is approximately periodic, as in ID simu- 
lations, but, in case 'a', we obtain recurrent bursts of time 
duration not much smaller than the cycle-time (of about 2 
s), whereas in ID simulations the bursts duration is very 
short in comparison with the cycle-time. 

Moreover, convection and circulation also affect the un- 
stable zone extension: this zone is not the full A zone, ex- 
tending up to 30-Rg, but only a portion of it (up to 13-Rg), 
since the convective motion stabilizes a large fraction of the 
radiation pressure dominated region. 

In case 'b', the low and high states have about the same 
duration (one half of the cycle-time, that is, on average, 50 
s). In ID simulations, instead, the high state has an ex- 
tremely short duration. 

Also, we have compared the disc structures we ob- 
tained from our simulations with the corresponding con- 
figurations calculated from the one-dimensional Shakura- 
Sunyaev model. The agreement between the results of the 
two approaches (the 2D simulated model and the ID theo- 
retical one) is in general good (see figs. 2 and 4 for case 'a' 
and 8 and 10 for case 'b'), taken into account that we have 
a time varying disc structure. The only significant differ- 
ences regard the innermost disc regions: for case 'a', between 
10 and 20 Rg, the thickness calculated from the canonical 
model is about 2 times the one of the simulated disc, and 
also the theoretical ID temperature is larger than the sim- 
ulated one (by 30 % about); for case 'b', from 3 to 15 Rg, 
the thickness calculated from the canonical model is about 
2-3 times the one of the simulated disc, while the theoretical 
ID temperature is larger than the simulated one by 30-35 
% about. So we confir m the result of Hure and Galliano 
jHure fc Galhandl200lll . who compared vertically averaged 
(i.e. ID) models of accretion discs with the corresponding 
vertically explicit (i.e. 2D) configurations and found that the 
one-dimensional description of the accretion discs structure 
is very close to the two-dimensional one. The generally good 
agreement between the disc structures we simulated and the 
Shakura-Sunyaev model allows us to conclude that the dif- 
ferent features between the thermal instability we found in 
2D and the limit-cycle behavior obtained by the ID calcu- 
lations are not due to differences in the used disc configura- 
tions. 

3) There are two time-scales that affect the time fea- 
tures of the limit-cycle phenomenon: the thermal time-scale, 
that determines the development rate of the thermal insta- 
bility, and the viscous time, that is connected to the A zone 
refilling after the collapse due to the instability. Here we 
discuss the role of these time-scales for the case 'a'. In this 
disc the collapsing zone extends nearly from 3i?g to 13-Rg. 
We give the values of the thermal and viscous time-scales 
at three points in this zone. These time-scales are quantities 
calculated from the theoretical formulae using the values of 
the disc physical variables resulting from the simulation. The 
expression we use for the visco us time-scale is t^isc = I v 
and, for the thermal time-scale dShakura fc SunvaevlllQTB) : 



itherm 



A{f3r] 



an 6{5l3r ~ 3) 
where A(l3r) = 8 + 51/3r 



(27) 



3/?,^ and f3r = Prad/P' 



Since the viscous time depends on the disc thickness, 
it assumes different values in the low (collapsed) state and 
the high (refilled) one. So we will distinguish its values in 
the two states using the labels 'LS' (Low State) and 'HS' 
(High State). 
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It is easy to see that the theoretical viscous time-scales 
are too large compared to the 'experimental' luminosity 
cycle-time, whereas the thermal time-scales are too small. 

We propose the following explanations. When the disc 
is in the low state, the accretion rate in the collapsed zone 
is small, but, at the outer boundary of this region, the disc 
is not collapsed and has a higher accretion rate. Because of 
this fact, the accretion flow at the outer boundary of the 
collapsed zone 'forces' matter to enter the region at a rate 
larger than the rate due to the viscous time-scale computed 
inside the region itself. So the refilling process is accelerated 
and its time-scale reduced with respect to the viscous time. 
A rough estimate of this effect can be given by the value of 
the viscous time-scale of the disk just outside the collapsed 
zone. Moreover, the calculation of this viscous time-scale has 
to take account of the 'real' viscosity present in the disc: the 
gas 2D motions (convection and circulation) give rise to a 
turbulent flow and consequently to an additional viscosity, 
of the order of magnitude Hvturb (where Vturb is the speed 
of the turbulent flow), that has to be added to the Shakura- 
Sunyaev one. 

To understand this procedure it must be remembered that, 
supposing a viscosity due to a turbulence (the Shakura- 
Sunyaev a- viscosity), we have not obtained a regular flow 
with the turbulence 'hidden' in the a-viscosity term. The 
flow produces another turbulence, not included in the a- 
viscosity term. Therefore, to express the total kinematic vis- 
cosity, we have to sum the term given by the speed of this 
new turbulence to the standard a-term. In formulae, we have 
a total kinematic viscosity given by: 



avsH + Hvtur 



(28) 



With this expression for the viscosity, the viscous time- 
scale tviac is given by: 

_ r2 _ r2 

tvisc — ~ ^JT- — (29) 

U aVsH + HVturb 

An approximate estimate for Vturb can be obtained by 
adding the gas radial and vertical speeds. These formulae 
give for t^isc the value of 12 s. Finally, it has to be con- 
sidered that the transition from the low to the high state 
is also accelerated by the process of the radial diffusion of 
radiation, whose typical time-scale is At — , where Al 
and T are the length and the optical depth of the region 
through which the radiation diffuses. Taking account of 
all these effects, the characteristic time of the transition 
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LS-HS is lowered from 12 s to about 5, value not far from 
the luminosity cycle-time we obtain in case 'a' (2 s). 
As regards the inverse process, i.e. the transition from the 
high to the low state due to the thermal instability, we 
form the hypothesis that the instability development time, 
that is esse ntially the thermal t ime-scale, is increased by 
convection dShakura et alJ Il978ll . Convection is naturally 
simulated by our 2D code, whereas the ID codes, obviously, 
cannot track the convective motions of the fluid masses 
and therefore do not include convection. This could be the 
reason for which, in ID simulations, the high state duration 
is very short: the thermal time-scale is not increased by 
convection and therefore the instability develops very 
rapidly, causing the collapse of the high state in the low 
one in a very short time. 

To support this hypothesis refer to fig. 6 for the radial 
speed in the two states low and high. It is clearly shown the 
extremely low radial speed of the collapsed zone together 
with the close 'active' zone exhibiting larger radial speeds. 
An oscillatory radial behavior of Vr is also evident. 
It is interesting that Szuszkiewicz and Miller also found 
similar radial speed profiles (Szuszkiewicz & Miller 2001), 
with significant oscillations in the behavior versus r. They 
claim that these oscillations are a numerical effect and, 
as a proof of that, show that, if an artificial diffusion is 
introduced, the oscillations disappear. We think, instead, 
that the radial speed oscillations are a real physical phe- 
nomenon, connected to the gas circulation in the disc, and 
that an artificial diffusion can obviously smooth away the 
oscillatory behavior, but this is only an artificial result due 
to a not physical ingredient. 

Finally we want to highlight the differences be- 
tween our results and Sz uszkiewicz and Miller's ones 
iSzuszkiewicz fc MilleilEoO J) about the hot gas bulge gen- 
erated by the thermal instability. In their work this bulge 
is just the medium by which the instability propagates 
through the disc. When the instability arises, a small 
region near the black hole becomes thicker and hotter. Its 
thickness and temperature become larger and larger while 
the radial extension also increases. The hot gas bulge that 
is formed through this process reaches a radial extension 
of 90Rg, then it begins to go down in the region near the 
black hole. The cooling wave that accompanies this process 
propagates gradually towards the internal boundary, until 
the whole gas bulge has returned to the disc equilibrium 
values of thickness and temperature. In our simulations the 
situation is different. The instability starts with the inner 
region collapse (therefore with the gas cooling) and not 
with the thickness and temperature increasing (the inverse 
processes). Both the collapse and the subsequent refilling 
are approximately simultaneous over all the unstable 
region: no heating and cooling waves propagate through 
the disc. When the unstabl e region is refilled, something 
similar to the hot bulge of (Szuszkiewicz fc Miller! (1200 iT) 
is formed. The involved region swells, but much less than 
in Szuszkiewicz & Miller's simulation. It reaches a vertical 
thickness not much larger than the equilibrium value and, 
as said, does not expand radially. 



6 CONCLUSIONS 

We put forward evidence with true 2D simulations that the 
limit-cycle behavior produced by the Shakura-Sunyaev in- 
stability is present in accretion discs having a radiation pres- 
sure dominated zone (A zone) . The time-scale of the instabil- 
ity and the shape of the light curve depend on the accretion 
rate. Lower accretion rates produce shorter time-scales of 
the oscillations. 

The 2D real motions play an important role to calcu- 
late the appropriate values of the oscillation frequencies. We 
obtain, for the subcritical regime (case 'a'), a frequency u of 
about 0.5 Hz and, in the supercritical case (the 'b'), f ~ 0.02 
Hz. In general the 2D time-scales are shorter (and the fre- 
quencies higher) than the ID ones. We attribute this result 
to the shorter viscous time-scale characteristic of the zone 
outside the collapsed region and to the role of the large con- 
vection present in the high state. 

These results may be relevant for the explanation of 
QPO emission in black hole candidates. More refined mod- 
els are necessary for a detailed interpretation of QPO ob- 
servational data, which, however, is not the purpose of this 
paper. Our aim is not to find an explanation of some sources 
behavior, but simply to see what happens to the standard 
disc structure during its time evolution. 
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